library(SDMTools)
species.data=read.csv('/home/jc112725/CARABIDS/CARABIDS/species_list_with_range1990.csv')
species=species.data$Species; species=as.vector(species)

year="2080"
scenario=c("sresb1","sresa1b","sresa2")

#set working directory where the raster files will be pulled from
work.dir = '/home/jc112725/CARABIDS/models/'; setwd(work.dir)
#List the rasters which have .asc.gz
#create empty object
richness_b1= NULL
richness_a1b= NULL
richness_a2= NULL
#make  aloop opening first asc file
for (spp in species)){ cat(spp,'\n')
sresb1=read.asc.gz(paste(work.dir,spp,"/summary/",scenario[1],".",year,".mean.realized.asc.gz", sep=""))
sresb1[which(sresb1>0)]=1
if (length (richness_b1)==0){richness_b1=sresb1} else {richness_b1=richness_b1+sresb1}
sresa1b=read.asc.gz(paste(work.dir,spp,"/summary/",scenario[2],".",year,".mean.realized.asc.gz", sep=""))
sresa1b[which(sresa1b>0)]=1
if (length (richness_a1b)==0){richness_a1b=sresa1b} else {richness_a1b=richness_a1b+sresa1b}
sresa2=read.asc.gz(paste(work.dir,spp,"/summary/",scenario[3],".",year,".mean.realized.asc.gz", sep=""))
sresa2[which(sresa2>0)]=1
if (length (richness_a2)==0){richness_a2=sresa2} else {richness_a2=richness_a2+sresa2}
}
#write.asc(richness_b1,paste("/home/jc112725/CARABIDS/richness_",scenario[1], sep=""))
#write.asc(richness_a1b,paste("/home/jc112725/CARABIDS/richness_",scenario[2], sep=""))
#write.asc(richness_a2,paste("/home/jc112725/CARABIDS/richness_",scenario[3], sep=""))
richness_b1= read.asc(paste("/home/jc112725/CARABIDS/richness_",scenario[1], '.asc',sep=""))
richness_a1b= read.asc(paste("/home/jc112725/CARABIDS/richness_",scenario[2], '.asc',sep=""))
richness_a2= read.asc(paste("/home/jc112725/CARABIDS/richness_",scenario[3], '.asc',sep=""))

current=read.asc('/home/jc112725/CARABIDS/richness.asc')
pnts=cbind(x=c(146.1,146.25,146.25,146.1), y=c(-15.7,-15.7,-16.5,-16.5))
Colormap= c("deepskyblue4",heat.colors(100)[seq(90,1)])

png ("/home/jc148322/2080_richness_scenarios.png",bg= "white",height=800, width=2000)
 
par(mfrow=c(1,4),cex=1,mar=c(2,2,2,0),xpd=T)

         image(current, main= "Present",cex.main= 2.5,cex.axis= 1.5,cex.lab= 1,col=Colormap,xlab="", ylab="",axes=F)
         axis(1);axis(2)
         image(richness_b1, main= "B1- 2080",cex.main= 2.5,cex.axis= 1,yaxt="n",xaxt="n", cex.lab= 1,col=Colormap,xlab="", ylab="",axes=F)
         axis(1)
         image(richness_a1b, main= "A1B- 2080",cex.main= 2.5,cex.axis= 1,yaxt="n",xaxt="n", cex.lab= 1,col=Colormap,xlab="", ylab="",axes=F)
         axis(1)
         image(richness_a2, main= "A2- 2080",cex.main= 2.5,cex.axis= 1,yaxt="n",xaxt="n", cex.lab= 1,col=Colormap,xlab="", ylab="",axes=F)
         axis(1)
          legend.gradient(pnts,cols=Colormap,limits=c(0,15), title='Species Richness', cex=1)

mtext("2080",side=3,at=.11,cex=1.8,line=0,outer=TRUE)

   dev.off()




